New constructions of domain decomposition 
methods for systems of PDEs 

V. Dolean* F. Natal] G. Rapin* 
1 fevrier 2008 



Resume 

We propose new domain decomposition methods for systems of partial 
differential equations in two and three dimensions. The algorithms are 
derived with the help of the Smith factorization of the operator. This 
could also be validated by numerical experiments. 



1 Introduction 

Neumann-Neumann BGLTV89 or FETI type algorithms are very popular 
domain decompositions methods. They are currently used for very large scale 
computations, see for example DDO03 and references therein. These methods 
are very well understood for symmetric definite positive scalar equations. For 
nonsymmetric problems and systems of equations many questions are still open. 
We propose in this note, a systematic construction of related algorithms for 
systems of partial differential equations (PDE) . The approach is based on the 
Smith factorization of the system of PDEs. First, we explain the derivation of 
the domain decomposition method for the Stokes system. Then the application 
to the Oseen system |NR| and the Euler equations |DN05| is briefly discussed. 

2 Analysis of the Stokes system via Smith fac- 
torization 

Let v > 0, we write the 2D Stokes equations as : 



^-Stokes 









-1 








X I I v I = I f v I . (1) 
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We first recall the definition of the Smith factorization of a matrix with poly- 
nomial entries and apply it to systems of PDEs : 

Theorem 2.1 Let n be an integer and A an invertible n x n matrix with poly- 
nomial entries with respect to the variable X : A = {&ij(X))-L<ij< n - 
Then, there exist matrices E, D and F with polynomial entries satisfying the 
following properties : det(E)=det(F )=1, D is a diagonal matrix and A = EDF . 

More details can be found in {W RL9 ,5j . We first take formally the Fourier trans- 
form of system Q with respect to y (dual variable is k). We keep the partial 
derivatives in x since in the sequel we shall consider a model probem where the 
interface between the subdomain is orthogonal to the x direction. We note 

/ -v{d xx - k 2 ) -0 X \ 

Astokes = -v(d xx - k 2 ) ~ik . (2) 

V d x ik / 

We can perform the Smith factorization of Astokes by considering it as a matrix 
with polynomials in d x entries. We have 

Astokes = EDF (3) 

where flu = D 2 2 = 1 and D 33 = — vA 2 with A = d xx — k 2 . One should note 
that a stream function formulation gives the same differential equation for the 
stream function. In the same way, the three-dimensional case can be charac- 
terized. In this case, the diagonal matrix D 3 £, is a four by four matrix whose 
entries are : D 3D U = D 3D 22 = 1, D 3D 33 = -uA and D 3D 44 = -vh 2 . This 
suggests that the derivation of a DDM for the bi-Laplacian is a key ingredient 
for a DDM for the Stokes system. 

3 A Domain Decomposition Method for the bi- 
Lplacian 

Let Q be an open subset of R 2 and T — ti n {x = 0} be a symmetry axis 
of Q. For simplicity, in the note we assume homogeneous Dirichlet conditions 
on the boundary d£l. The domain is decomposed into Oi = fl H {x < 0} and 
£1-2 = H {x > 0}. We consider the following algorithm : 

Starting with an initial guess satisfying = and A(w'j ) ) = Afa®) on T, the 



correction steps are expressed as follows for i = 1,2 : 

-vA 2 {w>l +1 ) = Omft, (4) 

d ni Atf +1 = -{d ni Aw k i + d n2 Aw k 2 )/2 (5) 

and d^w^ 1 = -{d ni w\ + d n2 w\)/2 on Y (6) 

followed by an update step : 

-vA 2 {w k ; +1 ) = g in n 4 (7) 

Aw k+1 = A< + (A< +1 +A^' +1 )/2 (8) 

andu;f +1 = w k + (w^ +1 + w^ +1 )/2 on T (9) 
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By symmetry arguments, we have converge in two steps to the solution of 
— vA 2 (w) = g in fl. 

4 The algorithm for the Stokes system 

Thanks to the Smith factorization (J3J, it is possible to translate the above 
algorithm for the bi-Laplacian operator into an algorithm for the Stokes system. 
It suffices to replace equations by the Stokes equations and in the inter- 

face conditions © , © , © , © w by the last component of F(u,v,p) T . In order 
to write the resulting algorithm in an intrinsic form, we introduce the stress 
a(u,p) on the interface for a velocity u = (u, v) and a pressure p. For any vector 
u its normal (resp. tangential) component on the interface is u n (resp. u T ). We 
denote a n and a T the normal and tangential parts of er, respectively. The new 
algorithm for the Stokes system for the same geometry as above reads : 
Starting with an initial guess satisfying = T2 and a® — — <j\ on T, 
the correction steps is expressed as follows for i = 1, 2 : 

in (10) 

+4nJ/ 2 ( n ) 
-(M^) +a T2 (S 2 fe ,p 2 fe ))/2 on r (12) 

followed by an update step : 

/ in Qi (13) 
< n + (^+ 1 +^+ 2 1 )/2onr (14) 

a„ 2 (^ +1 ,g 2 fc+1 ))/2onr (15) 

The boundary conditions in the correction step involve the normal velocity and 
the tangential stress whereas in the update step they involve the tangential ve- 
locity and the normal stress. In 3D, the algorithm has the same definition. By 
construction, it converges in two steps. In the iterative version of the Neumann- 
Neumann algorithm for the Stokes system |TP97| , P W02 , the boundary condi- 
tions of the correction step involve all the components of the stress whereas the 
update step involves all the components of the velocity. It can be shown that the 
convergence in two steps is then lost. More precisely, one obtains a convergence 
rate of 1/3 in the case n = R 2 , cf [NH] . 

5 Algorithm for other systems of PDEs 

The derivation of the algorithm for the Stokes system is based on the use 
of the Smith factorization and on the existence of superconvergent algorithms 
for scalar PDEs. The same procedure can be performed for the Oseen equa- 
tions [NR . In this case, the diagonal form of the operator is the product of 



Astokes[Wi ,q t ) 
anda T4 (^+\g? +1 ) 



AstofceHUj ,Pi ) = 
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a convection-diffusion operator and of a Laplacian operator. Using ( |AN97| V 
ATNVOO ) for the convection-diffusion, it is possible to derive an algorithm for 
the fourth order problem that converges in two steps. Translating this algorithm 
on the system, we obtain an algorithm converging in two steps for the Oseen 
system. The same work was done for the compressible Euler system |DN05| . In 
this case, the diagonal form of the operator is a product of a convective operator 
with a Helmholtz convective wave equation. The Smith factorization has also 
been used to design PML for the time-dependent compressible Euler equations, 
|Nat| . 

6 Preliminary Numerical results for the Stokes 
system 

The domain il — (—A, B) x (0, 1) is decomposed into two subdomains £!i = 
(-A,0) x (0, 1) and fl 2 = (0,5) x (0, 1). We compare the algorithm of §Hto 
the iterative version of the Neumann-Neumann algorithm. The stopping criteria 
is that the jumps of the normal derivative of the tangential component of the 
velocity are reduced by the factor 10~ 4 . In tablc2](lcft) A = B = 1, we see that 
both algorithms are not sensitive with respect to the mesh size. Of course, due 
to the discrete approximation we cannot expect the optimal convergence in two 
steps. But we only need one more step to achieve the error bound. We have also 
varied the width of the subdomains, (middle table). As expected the convergence 
of the Neumann-Neumann method deteriorates. For large aspect ratios, the 
method diverges (- in the table), since there exists an eigenvalue of the operator 
corresponding to the Richardson iteration with a modulus larger than 1. But still 
in this case convergence can be enforced by its use as a preconditioner in Krylov 
method as it is usually the case. Our new algorithm seems to be surprisingly 
robust with respect to the subdomain widths. For moderate variations we always 
need 3 iterations steps. If we choose very thin subdomains, for instance A = 1, 
B = 20, the stopping criterion is achieved in only 7 steps. In table (right), 
we have added a reaction term c > to the first two equations of the Stokes 
system. For instance c may be the inverse of the time step in a time-dependent 
computation. We see that the new algorithm is fairly stable. 

References 

[AN97] Y. Achdou and F. Nataf. A robin-robin preconditioner for an 
advection-diffusion problem. C. R. Acad. Sci. Paris, 325, Serie 
I :1211-1216, 1997. 

[ATNV00] Y. Achdou, P. Le Tallec, F. Nataf, and M. Vidrascu. A domain 
decomposition preconditioner for an advection-diffusion problem. 
Comput. Meth. Appl. Mech. Engrg, 184 T45-170, 2000. 

[BGLTV89] Jean-Frangois Bourgat, Roland Glowinski, Patrick Le Tallec, and 
Marina Vidrascu. Variational formulation and algorithm for trace 



4 



h 


new alg 


N-N 




B 


new alg 


N-N 


0.02 


3 


10 




1 


3 


11 


0.025 


3 


12 




2 


3 


12 


0.05 


3 


11 




3 


3 


11 


0.5 


3 


11 




5 


3 


15 


0.1 


3 


11 




10 


3 




0.2 


3 


10 




20 


7 





c 


new alg 


N-N 


0.001 


3 


11 


0.01 


3 


16 


0.1 


3 


19 


1 


3 


19 


10 


3 


16 


100 


3 


10 



Tab. 1 - Comparison between the new algorithm and the Neumann-Neumann 
algorithm (NN) : Iteration counts for different mesh sizes (left), aspect ratio 
(middle) and different reaction terms (right) 
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